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ABSTRACT 

In this paper I consider the nonlinear evolution of a rare density fluctuation 
in a random density field with Gaussian fluctuations, and I rigorously show that 
it follows the spherical collapse dynamics applied to its mean initial profile. This 
result is valid for any cosmological model and is independent of the shape of the 
power spectrum. In the early stages of the dynamics the density contrast of the 
fluctuation is seen to follow with a good accuracy the form 

«y = (i-(y L /i.5)- 1 - 5 -i, 

where 5l is the linearly extrapolated overdensity. 

I then investigate the validity domain of the rare event approximation in 
terms of the parameter v = 5l/<j giving the initial overdensity scaled by the rms 
fluctuation at the same mass scale, and find that it depends critically on the shape 
of the power spectrum. When the power law index n is lower than —1 the departure 
from the spherical collapse is expected to be small, at least in the early stages of 
the dynamics, and even for moderate values of v (|z/| > 2). However, for n > — 1 
and whatever the value of v, the dynamics seems to be dominated by the small- 
scale fluctuations and the subsequent evolution of the peak may not be necessarily 
correlated to its initial overdensity. I discuss the implications of these results for 
the nonlinear dynamics and the formation of astrophysical objects. 
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1. INTRODUCTION 

The gravitational instability scenario is the most widely accepted answer 
to the problem of the formation of the present structures of the universe. In 
such a view the astrophysical objects such as galaxies or clusters of galaxies are 
thought to be formed by the collapse of initial density fluctuations existing in the 
density field after the recombination. If it is well known that such perturbations are 
gravitationally unstable and are bound to grow under the influence of gravity, the 
complete dynamics of a density perturbation embedded in a fluctuating density field 
is not understood. Indeed it requires the resolution of nonlinear equations which has 
not been done yet. However the behavior of particular objects for which the initial 
conditions are well defined, such as a perturbation with an initial top-hat density 
profile, can be calculated up to their final collapse. These results have given birth to 
some theoretical derivations of the mass distribution function of the objects present 
in our universe with the so called Press and Schechter formalism (Press & Schechter 
1974). The main hypothesis on which this formalism is based is that any part of the 
universe will eventually collapse and form a virialized object according to its initial 
overdensity and within the time scale given by the spherical collapse model (e.g., 
Peebles 1980). It has also been suggested that it is more likely that the objects form 
at the location of the peaks, maxima of the density field (Kaiser 1984). It leads to 
the peak formalism developed by Bardeen et al. (1986). In any case, however, the 
time scale required for the formation of a virialized object is thought to follow the 
spherical collapse model. These formalisms are the only analytical models available 
to derive the epoch of formation of astrophysical objects such as quasars (Efstathiou 
& Rees 1988), or to derive the number density of objects like clusters (Peebles, Daly 
& Juszkiewicz 1989). Except with the use of numerical simulations, this is the only 
way to constrain the various cosmological models with such observations. It is thus 
of crucial interest to know whether the spherical collapse is a good approximation 
for the formation of isolated structures in a fluctuating density background. 

Peebles (1990) recently reconsidered this old problem arguing that the 
existence of small-scale fluctuations should slow down the dynamics of the collapse. 
In his view the non-radial motions should initiate a "previrialization" process. He 
indeed presents some numerical experiments in which such a trend appears with a 
large magnitude. These results may be at variance with what has been previously 
obtained by Efstathiou et al. (1988) and by Bond et al. (1990) when they checked 
the validity of the Press and Schechter formula, and also at variance with the 
results of Little, Weinberg & Park (1991) and Evrard & Crone (1992) in which 
they explicitly check the influence of a variable cut-off at small scale in the power 
spectrum, that fails to detect any "previrialization" effect. 

In this paper I consider this problem from an analytical point of view by 
analysing the early stages of the dynamics of a constrained Gaussian random 
field with perturbation theories. This paper, however, is not devoted to a new 
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perturbative theory but aims to consistently calculate some statistical quantities 
describing the dynamics of a given density fluctuation. I focus my interest on the 
evolution of the volume of such a rare fluctuation during its nonlinear evolution. 
The volume is simply defined as the volume occupied at any time by the particles 
that were in the initial density perturbation. Although it is certainly a naive picture 
for the formation of an object that likely involves accretion or ejection of matter, 
such an approximation is in the spirit of the current theories of the mass distribution 
function for which it is assumed that the matter content of a fluctuation is conserved 
during the collapse. In the discussion I will comment on the relevance of this 
hypothesis. 

In part 2, 1 compute the expectation value of the size of a density perturbation 
in the limit of rare events. It is shown that in such a limit it can be derived even in 
the fully nonlinear regime and it then follows the spherical collapse solution. Part 3 
is devoted to the discussion of the accuracy of this limit, that is, the first corrections 
expected from this limit case and their magnitudes. The case of a realistic model 
for the structure formation, i.e. a CDM model, is considered. A large fraction of 
the mathematical content of the resolution of this problem is given in appendices. 
The first one is devoted to usual results on the expectation values of non-Gaussian 
variables, and the second to the computation of two integrals involving the top- 
hat filter function. The third one is the most important. It contains the rigorous 
demonstration of the results presented in part 2. The last one is devoted to tedious 
calculations of interest for the discussion presented in part 3. 



2. THE NONLINEAR EVOLUTION OF RARE EVENTS 

2.1 The initial conditions 

I assume that the structures of the universe form from gravitational 
instabilities created by Gaussian fluctuations in a pressureless matter density field. 
The overdensity field 5(r) is then a random field so that its Fourier transform 
components 5^ defined by, 



are random complex variables with independent phases. The reality constraint of 
5(r) implies 



where ^(k) is the three dimensional Dirac distribution. The power spectrum, P(k), 
determines the full process of the structure formation. 




(1) 
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Throughout this paper I assume that within a spherical volume of radius Ro 
there is a given initial overdensity Si, 

1 



d 6 r 8(r). (3) 

Vo 



If the constraint (3) is fulfilled I will say that there is a "peak" in the given volume 
with the initial overdensity Si. This is a rather poor definition that makes sense 
at least when Si is large in a sense that will be given afterwards. One could 
have thought of stronger constraints for a given location to be called a peak, and, 
following Bardeen et al. (1986), required that this volume is at a maximum of the 
density field. Such a requirement will not change the qualitative results presented 
throughout this paper as it will be discussed at the end. Moreover the results 
presented in the following are valid for positive values of Si as well as negative 
values. So it will describe also the nonlinear evolution of rare voids. 

Once a peak is supposed to lie at a given location the resulting constrained 
random field is not only a function of P(k) but also a function of S{. This quantity 
is a Gaussian variable related to the S-^. Its variance is given by 

af = (6f) = * / d 3 r dV<5(r)5(r')> 

JVo 

For convenience I introduce the Fourier transform of the top-hat window function, 

— / d 3 re ikr = —- [sm{kRo)-kR cos{kRo)} = W TIi {kRo). (5) 
Jvo [kRo) 

The resulting variance for the initial overdensity is 



2 i d 3 k 



< = ] ^ W^(kR )P(k). (6) 

The purpose of this paper is to study the expected dynamics of the random 
field once the constraint (3) has been set. One obvious and direct thing to do is 
to consider the expectation values of the Fourier components of the random field. 
In the following I denote the expectation value of any random quantity when 
the constraint (3) is set. A priori any such quantity depends both on the power 
spectrum and on the value of Si. The variables Si and S-^ are both Gaussian and 
they are correlated together, so that 

P(k)W T}i (kRo) 

(2tt)3/2 at 
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The mean initial field exhibits a symmetric density perturbation, the profile 
of which is determined by the shape of the power spectrum. The expected local 
overdensity within a distance r from the center of the perturbation is 

<* « ,-»„ = /Jwwywr) k (8) 

and the local overdensity at a distance r from the center can be easily derived from 
this formula. The nonlinear evolution of a spherically symmetric perturbation with 
this mean density profile is also straightforward to compute: the trajectory of the 
matter at a distance r of the center depends only on the matter contained inside the 
radius r. The evolution of the profile of the object is shown in Fig. 1 as a function 
of time for various initial power spectra. There are a priori no reasons for the 
nonlinear evolution of these mean profiles to reproduce the mean nonlinear profiles. 
However, for a very dense perturbation one would expect that the dynamics of the 
collapse somehow follows this behavior. This is the central problem addressed in 
this paper. 



2.2 The equations of motion 

In the nonlinear evolution of the objects with the previous profiles, the size of 
the perturbation turns out to be the only quantity that is independent of the shape 
of the power spectrum. This is the reason why, also in the nonlinear constrained 
dynamics, I focus my interest on it. The size of the perturbation is assumed to be 
simply given by the volume occupied by the matter that was initially inside the 
radius Rq all along the nonlinear evolution. A Lagrangian calculation is the most 
adequate approach for the purpose of this calculation. In the initial unperturbed 
field, fluctuations are assumed to create initial random displacements that will be 
amplified by the gravitational dynamics. The unperturbed comoving positions 
are denoted q, and the exact comoving positions, x, are related to q through a 
displacement field q), 

x = q+#(t,q). (9) 



The volume V(t) subsequently occupied by the particles in the real space is 
related to the displacement field through space derivatives. Indeed it reads 



V(t) 



d d x 



V 



V 



<9x 



<9q 



d 3 q. 



(10) 



The jacobian of the transformation J(£, q) defined by 



J(t,q) = 



(11) 
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is then the quantity of interest for the calculation of the volume, 
determinant of the matrix the elements of which are 



dxj 



= Sij + 



d qj 



It is the 



(12) 



where 8^ is the Kronecker symbol. The motion equation of the point being at the 
position x leads straightforwardly with the Poisson equation to 

V x .x + 2-Vx.x = -4ttG(p(x) -p). (13) 

CL 

The density p(x) at the position x is given by the matter conservation constraint 

|J(x)|p(x) = p. (14) 

In the following I assume that J(x) is positive so that |J(x)| = J(x). This 
assumption is valid in the quasilinear regime since J(x) is then close to 1, but 
breaks after the first shell crossing where J(x) is zero. With this assumption we 
then have 



a 



a 



J(t, q) V x .x + 2-Vx .x = -47T Gp[l- J(t, q)] . 



(15) 



We assume that the time ti at which the initial fluctuations start to grow is 
arbitrarily close to the origin. The decaying modes of any kind are then subsequently 
erased so that the only purely growing modes remain at the present time. This is 
likely to be an unjustified approximation for the study of the virialization processes 
but we are only interested here in the collapse prior to the virialization. The time 
dependence of the growing mode in the linear approximation will be denoted D(t) . It 
is just proportional to the expansion factor in case of an Einstein-de Sitter universe 
but in general depends on the cosmological parameters (see Lahav et al. 1991 for 
a review). Another consequence of crucial interest is that the rotational modes are 
suppressed, since they are simply expected to be diluted with the expansion (e.g., 
Peebles 1980). I then assume that the displacement field is non-rotational, that is 



V x x tf(x) = 0. 



(16) 



The resolution of the dynamics in such an approach is entirely driven by the 
behavior of the displacement field, so that the motion equations have to be written 
only in term of \&(t, q) and its derivatives. The equation (11) then reads 



J(t, q) =1 + Vq + 



(Vq.*) 2 -Z>*A« 



1 

+ 6 



(V q .*) 3 - 3V q .* £ ¥..,*,., + 2 £ SijSj.k**,, 



t-3 



ijk 



(17) 
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where Vq is the gradient taken relatively to q and the summations are made over 
the three spatial components. The equation (15) involves a time derivative operator 
that I denote T(.) and that is defined by 

T(A) = A + 2-A. (18) 

(X 

Expressed in term of the displacement field only the equation (15) reads 
T (V q .*) + V q .¥ T (V q .¥) - ]T vJ^T (%,) + \ (V q .*) 2 T (V q .¥) - 

-^r(v q .*)X;[*<j*i,* - v q .*r(* 4J ) +X;**j%* r (*M) (19) 

= -47rGp[l- J(t,q)]. 

The equations (17) and (19) do not allow a complete resolution of the dynamics. 
However, the displacement field is fully determined as soon as the non-rotational 
constraint (16) has been set. The resolution of these equations would lead to 
a complete knowledge of the nonlinear dynamics (before shell crossings). I am 
interested here in a particular quantity, i.e. the volume occupied by a peak all 
along its nonlinear evolution. 

2.3 The expectation value of the volume 

The volume of the object is obviously a non-Gaussian random variable and a 
complicated function of the random variables 5^, which are correlated to the initial 
linear overdensity S^. In this part I am interested in the expectation value of the 
volume, {V{t)) s , for a random density field that follows the constraint (3). This 
quantity is a priori a function of Vq, of the power spectrum P(k), of the initial 
overdensity Si and of time. Once the time and the shape of the power spectrum 
have been specified this is a function of Si and Uj that contains the normalization 
of the initial fluctuations. Equivalently it is a function of the linearly extrapolated 
overdensity Sl, 

h = wt) s " (20) 

where D(t) is the time dependence of the growing mode, and the linearly 
extrapolated rms fluctuations at the same scale, 

«■ = (21) 



The regime that is investigated concerns the collapse of rare peaks or the expansion 
of rare troughs before the first shell crossing. As a result the linear overdensity of 
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the object, 5l, is at most of the order of a few units but the rms fluctuation at 
the same scale, gl, is small. The rare peak limit then implies that the expectation 
value of V(t) has to be taken for a vanishing value of gl and a fixed value of 6l. 



A.16) 



Following the appendix A, such expectation value is given by (Eqs. A. 13, 



/+oo 
dzQ v (t, cr L ,5 L + iz) exp 
-oo 



r z 2 - 


r+oo 


r z 2 i 


/ dz exp 

J — oo 




I 2*1. 


. 2 al\ 



with 



Gv(t, G L ,y) 



g (V(t)S p L 



2p 



p=0 



P 



(22) 
(23) 



since 5l is a Gaussian variable. The mean volume is then a complicated function 
of 5l and gl which I want to calculate in the limit gl — > 0. The leading term in 
the numerator of the ratio (22) is given by a saddle point approximation (around 
z 0). The relation (22) then reads 



(V(*)> ,«M*,<tl = o,<s l ). 



(24) 



I would like to stress that this expression is exact in the rare peak limit but involves 
all orders of the perturbation theory through the function Qv(t, gl = 0,5l) (but 
assuming no shell crossing and irrotationality). If the rare peak approximation were 
released, the expression (24) would no longer be valid. The part 3 is devoted to an 
analysis of the accuracy of this limit, its real meaning, and its validity domain. 

Each term involved in the expression of Qv(t, gl, 5l) has then to be calculated 
when gl = (as we will see each term is actually finite in spite of the denominators 
appearing in Eq. (23)). The quantities of interest (y(t)8 p L ) c have then to be 
calculated at the lowest order in gl (here Sl is not a fixed parameter but a Gaussian 
random variable). 

The random variable V(t) can be expanded relatively to Sy,, 

oo 

v(t) = J2v (p) (t) 

p=0 

where V ^ (t) contains exactly p factors of 5^ in its expression. The magnitude of 
Vb\t) is 

V {p) {t) ~V 8 P L . (25) 

The order is simply the co moving volume Vo; the order 1 corresponds to the linear 
approximation: = —Vq5l', and so on. The cumulant, (V(t)<5£) , appearing 

in the expression of Qy (see appendix A for a definition) connects p + 1 random 
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variables and then requires the product of at least 2p random Gaussian variables. 
As a result the leading contribution to this cumulant when the rms fluctuation is 
small comes from the order p of V(t). These arguments are the same that lead to 
the scaling behavior of the p order cumulants of the distribution function of the 
density at large scale (Fry 1986, Bernardeau 1992). The leading contribution at 
small ol is the exact result in the rare peak approximation, so that in such a case 
we get 

j p (t) ee (V(t)8 p L )jV = (V^(t)Sl)jV (26) 
and this quantity is of the order of 

3 P (t) ~ of. (27) 

The function Qv(t, o~l = 0, Sl) is then a finite number for finite values of Sl- 

At this stage the most difficult part is the derivation of the function Gv(t, 0, y) 
(denoted Qy{t,y) in the following). It cannot be done by a simple perturbative 
calculation since all orders of the dynamics are involved. As it is a rather 
complicated and long demonstration I present it only as an appendix (appendix C). 
It contains absolutely no other approximations than the ones presented previously 
and it gives the values of the whole series of j p (t) as written in the equation (26). 

The result of the appendix C is quite simple and shows that the function 
Qy (t, y) follows the differential equation, 

d< (R r ,lGv(t,y)/V„} 1/3 ) 

which is the motion equation describing 

a spherical collapse, where (^Ro [Gv(t, 2/)/K>] 1//3 ) is the size of the object (A is a 
possible cosmological constant). The behavior of Qv(t,y) for small values of y, 
(Gv(t, y) = Vq[1 — y + . . .]) ensures that the corresponding overdensity is simply y. 
The function Qv(t, Sl) that gives the expectation volume of a peak of overdensity 
Sl is then simply the volume occupied at the time t by an object having the same 
initial overdensity and initial comoving volume, according to the spherical collapse 
model. This result is independent of the cosmological parameters (O and A). As it 
can be noticed, the dependence on the power spectrum completely vanished and an 
individual collapse is expected to follow the same dynamics whatever the shape of 
the power spectrum. This result, however, is not robust, as it is discussed in part 3. 
It is quite comforting to see that the spherical collapse has something to do with 
the collapse of an object embedded in a fluctuating universe. 



2.4 The nonlinear evolution of the profile, of the density 
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The previous calculation deals with the expectation value of the volume V 
occupied by the particles that were in the initial volume Vo- One obviously could 
have considered the particles that were initially within the radius XRo (whatever 
A) and derive the expectation value of the volume occupied by these particles. The 
same kind of calculations can be done for this more general problem and the result is 
that they occupy the volume A 3 Vq Gv[t, {$(< XRo)) s ] which is the exact nonlinear 
evolution of the mean profile (this result is not a consequence of the previous one 
but its demonstration is nearly identical). The result is also valid for underdense 
areas so that the spherical model solution is expected to describe as well the collapse 
of rare peaks as the expansion of rare troughs. 

I propose a simple analytical expression for the evolution of the volume valid 
as long as the density contrast is not too large: 



(nt)) SL = vo 



i — - 

1.5 



r n 1.5 
OL 



(29a) 



where 5l is the linear extrapolation of the initial overdensity (eq. [20]). The 
dependence with the cosmological parameters is entirely contained in the time 
dependence of D(t). The form (29a) is exact for O — > 0, A = and is a good 
fit for any values of £1 and A, especially for the underdense regions (5l < 0). It 
fails, however, to reproduce with a good accuracy the position of the singularity 
(Fig. 2). 

The mean overdensity of the object, simply given by (Vo/V(t)) g — 1, can 
also be calculated in the rare event limit. It involves a similar generating function 
Gs(t,y) that can be calculated in the rare event limit (appendix C). The result in 
such a case is quite simple and reads, Gs(t, y) = Vo/Gv(t, y) — 1, so that the mean 
density is simply given by the inverse of the mean volume and one can use the 
relationship, 



1 - 
1.5 



r -t -1.5 

OL 



- 1. (296) 



As a result it is shown that in the rare event limit the density fluctuations 
behave as if they had the mean initial profile described in the section 2.1. In 
particular the density fluctuations are expected to be spherical (otherwise the 
timescale of the collapse would not be the one of the spherical collapse) and the 
outskirts of clusters for instance are well described by such an approach (while the 
inner part is subject to processes that are beyond the scope of these calculations). 
This result may not be surprising. It is simply that the rare event constraint of the 
random Gaussian field (3) is strong enough to determine the shape of the initial 
profile and its subsequent dynamics. What has been obtained here is a rigourous 
demonstration that in the rare event limit, in a field with Gaussian fluctuations, the 
dynamics of a constained field follows the spherical collapse up to the first singularity 
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(that is at 5l ~ 1-68 for an Einstein-de Sitter universe). The subsequent evolution 
cannot be studied with a perturbation theory (the reason is that the determinant 
J(x) is not necessarily positive and taking the absolute value cannot be obtained 
by means of a perturbative expansion). 

In practice however, the astrophysical objects may approach such a limit but 
never reach it exactly. This is the problem addressed in the third part of this work. 



3. THE ACCURACY OF THE RARE EVENT LIMIT 

The interest of the previous calculations is not only that we recover a known 
dynamics but that it is now viewed as a limit process. The accuracy of this limit 
can now be checked for each step of the calculation. It can be evaluated either 
by the departure of the exact mean collapse from the spherical collapse, or by 
the magnitude of the scattering around such a behavior. Firstly I examine the 
first corrections to the spherical collapse that appear when the rare event limit is 
released. The parameter that is supposed to be large is defined by 

Sl 

v = — 

The probability distribution of v is normal and rare events correspond to large 
values of z/, either positive or negative. In case of clusters for instance the values 
for v derived from the observations are of the order of 2 to 4 (Peebles et al. 1989). 
For quasars the values required for v may be greater. In the following I will then 
be concerned by the values of v of this order of magnitude or greater. 

3.1 The corrections to the mean collapse 

The mean value of the volume can be expanded with respect to ctl, 

0% = <<> + *! «> + •■• 

The term (V)f^ is the one that has been calculated previously in the rare event 

limit. The first correction is given by the term cr L (V) s (the order 1 in ff L is 
identically zero) . To discuss the results in term of v this correction can be written 
as 5j j /v 2 (y) g which due to this change of variable gives a correction starting 

(2,) 

at the quadratic order in Sl. The evaluation of (V) s has to be done by a 
perturbative calculation around the rare event limit, that is around ol = 0. During 
the calculations of part 2 the o~l = limit has been invoked two times, once to 
derive the expression of the mean value of V(t) as given by (22), then to calculate 
j p , in Eq. (26). It would be possible to calculate (22) exactly but the derivation 



(30) 



- 12 - 



of the parameters j p is more complicated and cannot be done in general. Actually 

it will not be possible to derive (V) s exactly but only the first two terms of its 
expansion with respect to 8l- These terms are partly given by the first correction 
in o\ of jo and j±. Indeed in general the value of j p can be seen as an expansion 
with respect to 5^., 



+ 



(31) 



The first term of this expansion is the one already calculated, and it has 
been shown to correspond to the spherical collapse. The corrective terms cannot 
be calculated as easily and the calculation has to be made by hand for each of 
the terms. For p = the situation is extremely simple since the corrections are 
all identically zero. It just means that the ensemble average of the volume simply 
follows the co moving volume. The first non-trivial term appears for p = 1. The 
leading term is simply the one coming from the linear approximation and the first 
correction is then 

m = -al + y^\t)5 L ) c + ... (32) 

The calculation of this corrective term is presented in the appendix D in the case 
of an Einstein-de Sitter universe. According to the relation (D.19) the expression 
of ji is then 



1 - 



^Wfa(kRo)P(k) J ^-P{k')f{k>/k) 



= ~ a L [l ~ a L Cexp] 
where the corrective coefficient C e xp is given by 

[ I 0^TI# R ^ P ^ I 0zP(k')f(k'/k) 



(33) 



Cexp — 



21 



70 Wfa(kR )P(k) 



(34) 



and /(r) is the function defined in (D.20). It is plotted in Fig. 3. The corrective 
term appears to be the product of one non-dimensional parameter, C e xp, and a 
time dependent factor which is simply the linear extrapolated value of the variance 
at the scale Rq for the present time. This allows one to derive the exact form of 
{V(t)) s up to the third order in 5l as a function of v. The results are presented 
for an Einstein-de Sitter universe. In such a case we have 



M*) = ^i, hit) 



——erf. 
189 1 



(35) 
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The resulting form for the expectation value of the volume when the corrections 
coming from the exact calculation of the integral (22) are included, is 



(V(t)) 5 =Vo 



1 ; -3 + 567/10 Cexp 



v 2 



21 x ' ' u 567 

(36) 

The dynamics of the spherical collapse is recovered in the limit v — > oo. The 
accuracy of this limit then depends on the values taken by Cexp- As can be seen 
from the equation (34) the value of this coefficient is a function of the shape of the 
power spectrum. In case P(k) has a power-law shape 

P(k) oc k n (37) 

the resulting values for Cexp are given in Fig. 4a. The most remarkable feature 
is that the coefficient diverges for n > — 1. This divergence comes from the large 
values of k and is due to the shape of the function /(r) for r — > oo (/(r) ~ 5/8r 2 
when t — > oo). The small scale fluctuations are then responsible for the infinite 
value of Cexp when n > — 1, and for such power spectra the rare event limit can 
never be reached, whatever the value of v. In order to check the importance of this 
divergence, I also considered the case of a power spectrum with a cut-off k c for large 
values of k. I then give in Fig. 4b the behavior of Cexp when the position of the 
cut-off varies compared to the inverse of the size of the fluctuation. The results are 
given for various values of n. It appears that for n = —1, which is the critical case, 
the divergence is only logarithmic so that a very large dynamical range is required 
for such an effect to be noticed. For n = the divergence is more rapid (as [Rok c ] 2 ) 
and such an effect could be seen in numerical simulations. 

When n < — 1 the corrective coefficient is small and the values of v of a few 
units are enough to ensure a behavior close to the spherical case. For n = —2 we 
have Cexp ~ 0.11. In Fig. 5, I present the expected density of the object (the mere 
inverse of the volume) as a function of the initial overdensity for v = 4 and v = 2. 
The curves correspond to 



I 



^ 246/341 - 567/341 C eX p 



v 2 



8l+. 



(V /V(t)) SL = 1+^(1+^-)^ 

(38) 

The relation (38) has not been obtained as a mere expansion of the inverse of the 
relation (36) but as an exact determination of {Vo/V(t)) s up to the third order. 
It happens that it is the same. The departures from the spherical collapse (at the 
same order) remain small. 

For the known astrophysical objects this is then of crucial interest to know 
the value of n at their mass scale. For galaxies and quasars a direct observation 
is not possible since the universe is now fully nonlinear at their mass scale, but 
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according to most of the theories based on the gravitational instability scenario, 
n is likely to be of the order of —2 or smaller. For the clusters the situation is 
more critical. Direct measurements of the mass fluctuations through galaxy counts 
give n ~ —1.3 according to Peacock (1991) for the APM angular survey, n ~ —1.4 
according to Fisher et al. (1993) for the IRAS galaxy survey. The use of the 
temperature distribution of the clusters allows one (with the Press and Schechter 
formalism) to do a direct estimation of the shape of P(k) for the matter field. It 
leads, according to Henry & Arnaud to n ~ — (l-7^ ; 3 5) and according to Oukbir & 
Blanchard (1992) to n ~ — 2 independently of the density of the universe. It then 
seems that for the clusters it is also fair to use the spherical collapse model. 

The case of the CDM model is considered in Fig. 6. The coefficient Cexp is 
computed with 

P{k) K (l + ak + pk^+^fk 2 ) 2 ' (39) 

a = 1.71/, (3 = 9.0/ 3 / 2 , 7 = 1.0/ 2 and / = 4.0 corresponding to h = 0.5, = 1. 
(Davis et al. 1985). The corrective coefficient remains essentially finite for the 
masses of interest, so that the spherical collapse model is expected to hold for such 
a shape of power spectra. 



3.2 Fluctuations around the mean behavior 



The purpose of this part is to derive the value of (F 2 (t)) 5 — (V(t)) g to 
evaluate the magnitude of the scattering around the mean collapse. 

First of all we can calculate this expression in the same limit, i.e. the rare 
event approximation. According to the appendix A, this function is given by 



/+oo 
dz [Q^(t, 5 L + iz) + Q v i(t, 5 L + iz)] exp 
-oo 



oo 
+oo 



2*1 



dz exp 



2*1 



(40) 



with 



p=0 



2P 



pi 



In the rare event limit, and for the same reason as before we have, 



(41) 



r=l 



The latter expression can be calculated exactly. The result given in the appendix 
C is quite simple and reads, 

i i 2 



(42) 
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This result is only valid in the rare event approximation: the function Qy that 
appears in this expression corresponds to the spherical collapse and the first 
corrections to (42) are of the order of a\. The calculation of the expression (40) 
has then to be made in a consistent way. We obtain 



(V 2 (t)) SL =G v (t,5 L )+aig v2 (t,5 L ) 

o\ \Gv(t,S L )^g v (t,8 L ) + 



l 2 N 



dy 



Gv(t,8 L ) 



+ 



(43) 



In the first term Qy does not only correspond to the spherical collapse but also 
contains corrections of the order of o\ when the approximation (24) is released. 
The last term corresponds to the first correction of (40) around the saddle point 
position (obtained by an expansion of Q v (t,y) around z = 0). The first neglected 
terms are of the order of a\. At the same order we have 



( V (t)) 2 5 L = Gv{t^ L )-alg v {t,8 L )^G v {t,8 L ) + 



(44) 



The equations (42-44) imply that all the terms that are independent of &L 
or proportional to a\ vanish, leading to a scattering of the form of, 



(v 2 (t)) Sr -(nt)) 



1/2 



a l Cscat(ML) Vo, 



(45) 



at the lowest order in a^. The coefficient C SC at(^^) * s a finite function, which 
in case of an Einstein-de Sitter universe is a function of Sl only. This result just 
means that there is no scattering when the rare event approximation is valid: the 
magnitude of the scattering is just equal to the magnitude of the error between the 
spherical collapse model and the exact mean behavior. 



To be more convincing I derive the value of C sca ^ for Sl 
examination of the various contributions leads to the expression, 



0. A careful 



1 



v 2 *l 



1/2 



(46) 



It is calculated in the appendix D and plotted in Fig. 7 as a function of n in 
case of a power-law spectrum (eq. [37]). Once again the coefficient becomes infinite 
when n > —1. For n m —2, C sca ^ ~ 0.3 so that the scattering around 8l = is of 
the order of 



(v 2 (t)) Sr ,-(nt)) 



1/2 



8 2 T 

0.3 V L 



(47) 
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In such a case the fluctuations around the spherical collapse are expected to be 
small. 

The fluctuations around the mean density can also be calculated. They are 
given by 

\m SL -(m)V 1/2 - s ' 



c f — 



(48) 



with 



c - 1 

^SCat T/2_-2 



\ (vM®)' > - 2(v^(t) (V«(t)) 2 ) + < (V 2 )(t)) 2 > 



l((vM(t)) 2 6l) 2 -(vW(t)8iy 



1/2 



(49) 



For n ~ — 2 we get Cg Ca ^ ~ 1- The situation obtained there is sketched in Fig. 8: 
the overdense objects are expected to follow the spherical collapse model with a 
possible departure and a possible scattering of the same order of magnitude. Then 
the accuracy of the spherical collapse model depends on the value of v. It is not 
possible to determine the value of 5 up to which this approximation is valid since 
the corrective term is known only up to the third order in 8l- At this order the 
discrepancy between the dynamics of the spherical collapse and the real dynamics 
is seen to be weak as long as v is of the order of a few units (say 2, 3) (eqs. [36, 
38]). It implies that the spherical collapse model is relevant at least up to 5 ~ 4 
and for values of v of a few units. Whether this approximation breaks down when 
6 is larger than 4 is a problem not solved by these analytic calculations. 

When n > — 1, however, the situation is quite different, and the fluctuations 
are unbounded. It seems that the initial overdensity 5l does not determine at all 
the subsequent evolution of the object which rather depends on the small-scale 
fluctuations. The origin of this behavior is discussed in the next section. 



3.3 Discussion 

I want to stress that the results presented here take fully into account the 
application of a filter function on the density field at any order of the expansion. It 
represents a considerable improvement compared to previous results. For instance 
Kofman (1991) proposes a density distribution function based in the Zel'dovich 
(1970) approximation, thought to be valid when the fluctuations are still small 
that, however, does not take into account the filtering of the evolved density field. 
Bernardeau (1992) also presents the expected shape of the density distribution in 
a quasi-Gaussian approximation (that takes into account the nonlinear behavior of 
the density field) with the same kind of approximation. Both of these calculations 
were made for a density in a vanishing volume and appear as approximations to 
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the exact solution for a finite and large volume. In the calculation presented here, 
however, the filtering process has been fully taken into account. It turns out that 
the expected nonlinear evolution of a rare event follows exactly what would have 
been expected from the spherical model. 

It has been stated that the geometrical dependences contained in the second 
or third order terms may induce a slowing of the dynamics by tidal or shear effects 
(Juszkiewicz & Bouchet 1992) that would produce a sort of previrialization. The 
calculations I presented here take into account these geometrical dependences but 
contain no indications of previrialization, at least in the rare event limit (and as long 
as n < —1). The idea that the small-scale fluctuations may alter the collapse rate 
of an object was expressed by Peebles (1990). He presented numerical experiments 
where indeed the growing rate of the fluctuation was slightly lower compared to 
the spherical collapse model even for objects of final density of the order of 3. This 
result is at variance with what has been obtained here in the case of moderately rare 
events. Previrialization however may be present in the final stage of the collapse 
when the radius of the object starts to decrease for a not extremely rare event. 

The extent of the central problem solved in this paper, the expected behavior 
of a constrained random field, also deserves some comments. In the equation (3) a 
"peak" has been defined as a mere overdense part of the random density field. The 
result of part 2, however, holds even for more constraining initial conditions. One 
could have imposed that the volume is at an extremum of the density field (in case 
of a rare overdense volume, it is likely to be a maximum) . In such a case the initial 
random field can be seen as a non-homogeneous random Gaussian field for which 
the relation (1) has been replaced by 

/A A \ X<V+V\P<h\ k!.k 2 P(h)P(k 2 ) W(h)W(k 2 ) 

\ 5 k/k 2 / = Wi + k 2 )P(/ci) — 737; (50) 

f W*(k) P(k) 

where W(k) is the shape of the window function in the Fourier space used to 
determine whether a point is the peak or not (this is not necessarily the top-hat 
window function). The random field is still Gaussian and isotropic around q = 
so that the whole demonstration still holds. The corrections to the rare event limit 
and in particular the values of the coefficients Cexp and Cscat' however, are likely 
to change (and to be reduced). One could have also required that a greater volume, 
Vi, centered on Vo contains a certain overdensity, 5i in order to study the formation 
of objects in a rich, or a poor, environment as did Bower (1991) in the frame of 
the Press and Schechter formalism. The conditional random field is still Gaussian 
and isotropic so that the demonstration still holds, and in the rare peak limit the 
dynamics follows the spherical collapse model. The only change is that the mean 
value of 5l is not zero anymore but 

{Sl) = j'd 3 kw- H (*ij )fW * (51) 
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(with an unchanged variance), so that the distribution of 5l is shifted towards higher 
or lower values according to the sign of 81 . This is also of great interest for the Press 
and Schechter formalism, since as stressed by Bond et al. (1990) and Blanchard, 
Valls-Gabaud & Mamon (1992), the mass distribution function is based on an 
assumption on the dynamics of a somehow isolated density perturbation: the initial 
constraint is not only that a given overdensity is reached in a given volume but also 
that the surroundings do not form a higher overdense region. This local constraint 
may be especially important for the low values of v. This paper does not address, 
by far, the problem of the exact nonlinear dynamics with such initial conditions 
when v is small so it does not provide a justification for the Press and Schechter 
formula in its whole generality. The large mass behavior, however, corresponding 
to the high v limit received a strong dynamical justification. The demonstration 
is unfortunately not complete, since the virialization process remains beyond the 
scope of these calculations. The work of Thomas & Couchman (1992) suggests that 
the spherical collapse model is correct up to the maximum expansion of the object 
but fails to reproduce the contraction rate. This is not that surprising, since in this 
regime the decaying and rotational modes and the shell crossings that have been 
neglected are likely to play a major role in the process of relaxation. 

The second major phenomenon highlighted in this paper is the existence of 
a critical value for the power-law index. The existence of an unbound correction 
when n > — 1 is the signal that a perturbation approach is not safe. Such a result 
deserves a deeper examination. It turns out that the divergence of C e xp is entirely 
due to dynamical properties. For instance if we consider the conditional expectation 
values of the spatial derivatives of the displacement field, we get, 



which is finite only if P(k) <C k~ x when k — > oo (since /(r) ~ 1/t 2 for r — > oo). 
As a result when n > — 1 the displacement field is not a derivable quantity: two 
points that are initially very close may have completely different trajectories. This 
is the indication that the matter follows a chaotic behavior rather than a mere 
coherent flow towards the overdense regions. It is confirmed by the derivation of 
the expectation value of the functional derivative of the displacement field with the 
initial fluctuation. The result is 




(52) 



*tti(q) v_ f 

5 5 k t J 



(2tt) 3 / 2 k 2 



d^k -ik 



1 e iq.k 



(53) 
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which is also divergent for n > —1. As can be seen in (53) this is only due to the 
shape of /(r) and not at all to the shape of the filter function. As a result a small 
change in the initial fluctuations would induce a dramatic change of the behavior of 
the displacement field. It seems that the particles that were initially in the volume 
Vq do, at least in the first stage of the dynamics, a sort of chaotic diffusion, and 
have no reason to form ultimately a unique virialized object. This effect is present 
whatever the value of n but when n < — 1 the correlation length of the displacement 
field becomes infinite (as it is for the velocity field) and the displacement affects 
the whole object without changing its volume. When n > — 1 this is no longer the 
case and small-scale displacements can disrupt the collapsing peak. That does not 
mean that no object forms at the position of the peak, but it is a complete nonsense 
to try to derive its geometrical properties with the trajectories of the particles 
that were initially inside Vq. The edge of the volume occupied by the particles 
is expected to be extremely fuzzy (like a fractal set), and the approach adopted 
here is totally inadequate to deal with this problem. The use of more constrained 
initial conditions (such as in eq. [49]) is not expected to fix this problem, since 
these constraints do not involve the small-scale fluctuations. A systematic study 
of the dynamics of constrained fields with numerical simulations would be helpful 
to clarify the situation, but in any case the simple picture used to derive the mass 
distribution functions, based on the assumption of mass conservation, is severely 
altered in the n > — 1 case. 
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APPENDIX A: Conditional expectation values of non-Gaussian variables 



This appendix is devoted to a recall of the general form of the distribution 
functions of Gaussian and non- Gaussian random variables once the series of their 
moments is known. We consider two random positive variables V\ and V 2 for 
which their high-order correlations and cross-correlations, (VfV^) p,q integers, 
are known. These moments define a two variable generating function 

M(\uX 2 )= £ (VfU)^. (Al) 

p=0,g=0 P ' q ' 

Then the two variable distribution function of V\ and V 2 reads 

/+ioo+0 i\ p+loo+0 11 
/ ^ 
-ioo+0 2nl i-ioo+O 2nl 

x M(Xi, A 2 ) exp(-AiVi - A 2 V 2 ) 



(A.2) 



Such a form can be checked by deriving the moments of the distribution 
function defined in (A.2). Indeed the expectation value of VfV^ is given by 



= (v?v,]). 



-+1OO+0 



QP QQ 



d\ p x d\ q 2 



Ai=0,A 2 =0 



(A3) 



However it turns out in fact that it is more efficient to work with the 
cumulants of the distribution, (VfV^)^ They are defined recursively from the 
moments: 

<i> c = °; 
(v) c = (v); 

<ViV 2 ) c = <V,K 2 > - (V 1 ) c <V 2 ) c ; 
{V 3 ) c = {V)-(V)l-3(V) c {V\; 



(A4) 



In general the cumulant of order (p, q) is obtained in such a way: consider all 
the partitions of a set of p + q points (except the one constituted by the set 
itself), take the product of the cumulants of each of the subsets that come from 
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moments of lower order, then the cumulant is obtained by substracting all the terms 
obtained such a way from the value of the moment of order (p, q). Note that the 
usual arithmetic operations are not allowed on the cumulants and for instance that 
(ABC) c 7^ ( (AB)C) c since the first is the cumulant of three random variables 
and the second of two. The cumulants also define a generating function %(Ai, A2) 



OO ,p > q 



p=0,q=0 



The functions %(Ai, A 2 ) and -M(Ai, A 2 ) are closely related together by the relation 
(e.g., Balian & Schaeffer 1989), 

^(A 1 ,A 2 )=exp[ x (A 1 ,A 2 )], (A6) 

so that the distribution reads now 

/+ioo+0 i \ r+loo + n 

-4 / T ^exp[ x (A 1 ,A 2 )-A 1 y 1 -A 2 V 2 ]. 
-ioo+0 ^ 7n </-ioo+0 ^ 7ri 

(A7) 

The general form (A. 7) allows one to calculate any statistical quantity related 
to the two variables V± and V 2 . One quantity of interest in this paper is the mean 
value of one knowing the other, (Vi)^ . This quantity is given by the ratio, 



/» 00 poo 

( Vl )v 2 = J dV i y i v(Vi,V 2 )/ J dV 1 V (Vi,V 2 ) 
" +ioo+0 dA 2 d 



o ■ x -x(Ai, A s 

ioo+0 ^7T1 CAl Ai=0 



exp[ x (0, A 2 ) - A 2 l/ 2 ] (A.8) 



/ / -^exp[ x (0,A 2 )-A 2 y 2 ]. 



-+100+0 

-ioo+0 

The partial derivative that appears in the result is given by 

a 



<9Ai 



x(Ai, A 2 ; 



We then apply the relation (A. 9) when V 2 is a Gaussian variable of 
expectation value 1, and of rms fluctuation a. This is the case for the early 
cosmological density field, so that we define the random variable 5l = V 2 — 1 for 
such a case. The function x(0, A 2 ) then takes a simple form since 



Xg(0,A 2 ) =X 2 + -a 2 X 2 2 



(A.10) 
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(the lowscript G is for Gaussian) . We can also notice that 

<F 1 (1 + ^> C = (^> C . 



(All) 



The expectation value of V\ knowing 5l is then given by the relation (A. 8) and it 
reads in this particular case (we make the change of variable y = a 2 A2), 



l +loo+0 dy r 



loo+O 



(y 2 - 2 yS L ) 



with 



GvM = Y, 

p=0 



2 a 2 



J-ioo+o 27ri 



exp 



(y 2 -2yd L ) l 
(2 a 2 ) \ 
(A12) 



I 



a 2p p\ 



(A.13) 



We can also compute the expectation value of V 2 knowing Sl- It is given by 
the following expression, 



,+100+0 d 

J -loo+O Z7V1 1 





-(y 2 - 2 yd L )] 


exp 


2 a 2 



+100+0 



lcx)+0 



dy_ 
2yri 



7 exp 



(y 2 - 2 yS L ) 



(AAA) 



with 



p=0 



(2 ° 2 ) 



a 2p p J 



(A15) 



The calculation of the integrals (A. 12) and (A. 14) can be simplified by the 
change of variable, y = 5l + iz. As a result we have 



/+00 
dz g Vl (5 L + iz) exp 
-00 



r z 2 ' 


r+00 


\ ^ 1 


/ dz exp 
J —00 


2a 2 _ 


2 a 2 . 



(A16) 



and 



dz 



/+00 
dz exp 
-00 



GI^Sl+iz) +g V 2(5 L + iz) 

~2 





r z 2 1 


exp 


2a 2 _ 



2 cr 



(A17) 
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APPENDIX B: Geometrical properties of the top— hat window function 



The purpose is this appendix is to give general properties of the top-hat 
window function, Wpj^x). It is defined by 



W TH (x) = 1 if |x| < R ; 
Wpjj(x) = otherwise; 

for a scale Rq. The Fourier transform of this function is then given by 



(B.l) 



W TR (kR ) = 



(sm(kR ) — kR cos(kR )) 



(B.2) 



where k is the norm of k. The latter function will be widely used for the calculation 
of the dynamics. 

Consider three wave vectors ki, k 2 and k 3 on which you may have to 
integrate. The properties of interest concern the integration over their angular 
parts, dfii, dfi 2 and dfi 3 , and are the following: 



and 



/ 



dOid0 2 W TH (|ki +k 2 | Rq) 



1 - 



ki-k 2 

k\k 2 



(4tt) 2 - W TH (k 1 Ro)W TR (k 2 R ) 



(B.3) 



dO!d0 2 d0 3 W TR (|ki + k 2 + k 3 | R ) 



x 



ki.k 2 

kik 2 



k 2 .k 3 

k 2 k3 



k 3 .ki 

k 3 ki 



+ 2- 



ki .k 2 k 2 .k 3 ka.k 



h 2 h 2 h 2 
n 'l n ' 2 n '3 



= (4tt) 3 - W TR (k 1 Ro)W TR (k 2 Ro)W TR (k 3 Ro). 



(B.4) 



These properties are essentially some properties of factorizability of the top- 
hat window function and are only true for this particular window function. It means 
that the top-hat window function commutes with certain geometrical quantities and 
that we can replace the angular dependence of k 1? k 2 and k 3 by their mean as if there 
was no window function at all. For instance the mean value of 1 — (ki.k 2 /fci/c 2 ) 2 
in (B.3) is 2/3 in the absence of window function which exactly corresponds to the 
coefficient 2/3 appearing in (B.3). 
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The proofs of these properties rely on the fact that WpjjW can be written 
in term of a Bessel function: 

W TR (k) = 3^/2 k~ 3 / 2 J 3/2 (k) (5.5) 

(the value of Rq has been set to 1 for simplicity but the demonstrations obviously 
hold for any value of Rq). The window function can then be written (e.g., 
Gradshteyn & Ryzhik, Eq. 8.532.1), 

°° 3 d 
W TH (|k 1 +k 2 |)=37r^(- + m)(A; 1 A; 2 )- 3 / 2 J 3/2+m (ki) J 3/2+m {k 2 )—P m+l {-u) 

771 = 

(B.6) 

where u = ki.k 2 //ci/c 2 and P m is a Legendre polynomial. 

Then I\{k\, k 2 ) defined as the left part of the relation (B.3) reads 



oo 



h(k u k 2 ) =3nJ2d + ™)(kik 2 )- 3/2 Js/2+m(ki) Js/ 2+m (k 2 ) (4tt) (2tt) 



2 

777 = 



X 



r 1 d 

/ du (1-u 2 ) —P m+1 (- U ) 

1 (5.7) 



3 

= -3tt ^(- + m)(/ci/c 2 )" 3/2 ^3/2+777(^1) ^3/2+777(^2) (4tt) 



777=0 

x y dtiti P m+ i(-u). 

The end of the demonstration is straightforward since u = Pi(u), and the Legendre 
polynomials are orthogonal to each other so that m = is the only non vanishing 
case. The latter gives j_ x duPi(u) =2/3, so that the integral reads 

hih, k 2 ) = (4tt) 2 3tt (hh)- 3 / 2 Js/^h) J 3/2 (k 2 ) (5.8) 

which, using the equation (B.5), gives the relation (B.3). 

The proof of the relation (B.4) is based on the same kind of calculation but 
the demonstration is more complicated due to the fact that three different vectors 
are involved. Two steps are required. First of all let me define K = ki + k 2 . The 
first step will be to separate K from k 3 . I can then define the angles 6i 2 , ©1, 2 
and O3 which are respectively the angles between ki and k 2 , between K and ki, 
between K and k 2 and between K and k3 (note that Q± 2 = ©1 + 2 ). The vector 
k3 is also defined by its recession angle, $3, around the vector K (they are set to 
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for the vectors ki and k 2 ). The angular part appearing in the left part of (B.4) is 



A(ei,e 2 ,e 3 ,$ 3 ) = 



1 - 



ki.k 2 



k 2 -k 3 

k 2 k? 



+ 2 



kik 2 J \ ^2^3 
kj.k 2 k 2 .k 3 k 3 .ki 



k^ki 

K /c 3 /ci 



^1^2 ^3 

= i-cos 2 (e 12 ) 

- [cos(e 3 ) 003(00 + sin(e 3 ) sin(6i) sin($ 3 )] 2 

- [cos(e 3 ) cos(6 2 ) + sin(e 3 ) sin(e 2 ) sin($ 3 )] 2 

+ 2 cos(0i 2 ) [cos(6 3 ) cos(6i) + sin(6 3 ) sin(6i) sin($ 3 )] 
x [cos(0 3 ) cos(0 2 ) + sin(6 3 ) sin(6 2 ) sin($ 3 )] . 

(5.9) 

The integration over $ 3 can be made at this stage since the factor 
Wpjj(|ki + k 2 + k 3 |) is independent of $3. The result reads 

s(0i, e 2 , e 3 ) = / d$ 3 a{q u e 2 , e 3 , $ 3 ) = 

Jo 



2tt ( 1 - cos 2 (6i 2 ) - cos 2 (6 3 ) cos 2 (6i) - - sin^(6 3 ) sin 2 (6i) 
cos 2 (e 3 ) cos 2 (e 2 ) - X - sin 2 (e 3 ) sin 2 (e 2 ) 

1 



cos 2 (6 3 ) cos(0i) cos(e 2 ) + - sin 2 (6 3 ) sin(6i) sin(6 2 ) 

2 



+2cos(6i 2 ) 

The function I 2 (/ci, /c 2 , /c 3 ), defined as the left part of (B.4), then reads 



(5.10) 



>{kiMM = J dfii d0 2 sin(e 3 )de 3 



x Stt £(- +m)(Kk 3 )-V 2 J 3/ 2+m(K) J 3/ 2 +m (ks) (5.11) 



m=0 



x — p m+ i(-cos(e 3 )) 5(61,62,63). 



The integration over 63 takes advantage of the fact that 5(0i,6 2 ,63) contains 
only terms independent of 63 or proportional to cos 2 (63). We can then use general 
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properties of the Legendre polynomials, 



f 1 d 

/ du— P m+ i(— it) = 2 for even values of m; 
J-i dw 

/_ 



d 

d-u— P m+ i(— it) k 2 = 2 for even values of m and m^O; (-^-12) 
! du 



/•l A 

/ dw— Pi(-u) w 2 = 2/3; 
7-1 d-u 



and these integrals are zero otherwise, for u = 003(63) to compute (B.ll). When 
m 7^ 0, the two contributions coming from the terms independent of O3 and the 
ones proportional to cos 2 (©3) have to be added which is obtained by setting ©3 = 
in S(Oi, ©2, ©3) which gives zero. 

For m = the contribution of the integral has to be calculated directly. 
Taking advantage of the fact that ©12 = ©1 + ©2 we obtain 

J sin(© 3 )d© 3 -^Pi(-cos(© 3 )) £?(©!, © 2 , ©3) = y (l-cos 2 (© 12 )) 

so that 

/Q A 
df2i d0 2 3n-(Kk 3 )- 3 / 2 J 3/2 (K) J 3/2 (k 3 ) — [l - cos 2 (© 12 )] 

= W T -fi(k 3 ) h(k 1 ,k 2 ) 

(B.13) 

which using the equation (B.3) gives straightforwardly the relation (B.4). 
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APPENDIX C: Calculation of Gv(t,y), Gv*{t,y) and G 5 {t,y) 

The function Qv{t,y) is denned in equation (23) and will be calculated for 
<jl = that is with the approximation (26) for the series of the coefficient j p . Let 
me first sketch the principle of this calculation. The starting point is the dynamical 
equations in Lagrangian space describing the evolution of a pressureless fluid (eqs. 
[C.2,C.3]). Then they will be written in Fourier space displaying the geometrical 
dependences examined in the previous appendix. Taking the cumulants of each 
term (once multiplied by Sl at the required power) introduces integrations over the 
angles of the wave vectors that are then known from the result of the appendix B. 
As the geometrical dependences vanish we end up with a system of equations of 
which G(t,y) is solution. The field equation system have thus been replaced by a 
simple first order differential system. 

C.l Notations and general properties 

The fulfillment of this derivation requires the definitions of new quantities 
presented in the following. Let me however start with a general property of the 
cumulants widely used throughout this appendix. In this paragraph any quantity 
denoted stands for the order p of Q relatively to the random variables 5^. The 
property of interest concerns the cumulant of a product. If Q can be written as a 



product, Q = F l F 2 , the order p of Q is given by = YX=o F i F 2 and the 
cumulant it gives when multiplied by 8 P L is 



This result is due to the properties of the Gaussian variables. A product of 2p 
Gaussian variables has an expectation value that is written as a summation of 
p moments of pairs, so that the terms involved in (C.l) separate in two sets: the 
factors involving F\ and the factors involving F 2 . The symmetry factor p\/r\/(p— r) ! 
is due to the number of possibilities you have to separate a set of p points into two 
sets of r and p — r points. 

The starting point of the derivation of Gv(t, y) is the motion equations derived 
in the main section, eqs. (17, 20), that I recall here 




(C.l) 




(C.2) 



+ * (Vq -*) 3 - 3V q vI '<,/ vI ',.< + 2 E 



ij ijk 
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where Vq is the gradient taken relatively to q. The second equation is 

T (V q .*) + Vq T (V q .¥) - £ tf^T + ± (V q 2 T (V q .¥) 

-^T (Vq £ - Vq T (*,,,) J + J2 *i,J*i,kT 

= -47rGp[l-J(*,q)] 

(C.3) 

where T(.) is a time derivative operator, 

T(A) = A + 2-A (C.4) 

The displacement field has also to be irrotational in the real space. These equations 
completely define the time evolution of the displacement field. A calculation of 
the first orders when the displacement field is considered as a small quantity is 
given in the appendix D. But such direct perturbative calculations are inadequate 
to compute the values of j p (t) as defined in (24). The series of the parameters will 
be obtained as a functional equation for Qv(t, y)- 

For convenience I define few other quantities. Similarly to Qy I can consider 
the function 

y^^l^W. (a5) 



P =0 *L P [ 



and the parameters 



d P (t) = ((±J d 3 qVq.v|>(t,q)^> c = ( (± J d 3 q Vq .^(t))^, (C.6) 

in the rare peak approximation. I also denote (t) and d^\t) the Fourier 
transform of respectively J(t, q) and V q .\I/(£, q) at order p in an expansion with 
respect to the random variables 5^. All these quantities are time dependent and 
the operator T(.) defined in (C.3) can apply to any of them. 

The term of order p for the volume then reads 

y(P)(t)/F ° = / (2^ W TR(kRo) j£\t) (C.7) 

and we have 

j > = I (2^^) ^> c ^TH(^o). (C.S) 
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We have similarly for the divergence of the displacement field, 



1 3i 



and 



^=/^< d k ) (*)^>c W TH(^)- 
The components of the matrix j also define some cumulants, 

Sp;iJ (t,q) = (*g(t,q)^> c . 



(C.10) 



(an) 



We can notice that the parameters j p (t) can simply be expressed with the quantities 
Sp;i,j since 

m = ^F [ d 3 q^ %l ,,(t,q). (C.12) 



Written in k-space the equations (C.2) and (C.3) introduce the functions 



Ji(ki,k 2 ) = W XH (|ki + k 2 | Rq) 



ki-k; 

kik 2 



(C.13) 



and 



J 2 (ki,k 2 ,k 3 ) = W TH (|ki + k 2 + k 3 | Ro) 



x 



1 - 



ki.k s 

kik 2 



k 2 .k 3 



k 3 .ki V 2 ki.k 2 k 2 .k 3 k 3 .ki 

/c 3 /ci 



h 2 h 2 h 2 



These functions have been studied in the appendix B where it is shown that 

2 



/ 



dQx d0 2 Ji(k!,k 2 ) = (4tt) 2 -W XH (A;ii2o) W TH (A; 2 i?o) 



(C.14) 
(C.15) 



and 



y dOi d0 2 d0 3 J 2 (k l5 k 2 ,k 3 ) = (An) 3 ^W TH (A;ii2o) W TR (k 2 R ) W TR (k 3 R ). 

(C.16) 

C.2 The function Qy(t, y) 

The first step of the demonstration is to remark that the quantities s p -ij(t, q) 
are symmetric in i, j and are functions of |q| only. The fact that the initial conditions 
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are randomly isotropic around q = implies directly that s p -ij is a function of the 
norm of q. Moreover when i ^ j, s p -ij — Sp-j^ are the components of a vector (given 
by the rotational of the displacement field) so that these quantities should be zero 
for the same reason: no particular direction exists in the initial conditions. Note 
that this property is correct at the level of these expectation values but not at all 
for the random displacement field itself. The random variable — ^ j^ is not zero 
although some of its expectation values vanish. The implication of such property is 
that s p -ij can be written 



We can now transform the equations (C.2) and (C.3) involving the 
displacement field in two equations between the series j p and d p . 

Let me start with the equation (C.2). This equation has to be true at any 
order so that for p ^ 0, 




(C.17) 



and, according to (C.12), s(k) is then simply given by 



s p (k) = (d^(t)Sl) c 



(C.18) 





r=p,s=p—r r- 



(C.19) 



r=0,s=0 L 



^ v q • i,j 3,i ^ L^i i,3 3,k k,i 



ij ijk 



The equation (C.19) can now be transformed in an equation between the 
cumulants j p and s p; ij by multiplying it by 8 V L and taking the cumulant of each 
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term. Using the equation (C.l), (C.17) and (C.18) we obtain 



ik.q 



1 p 



pi 



d 3 ki d 3 k 2 



2 ^ r\{p-r)\ J (2tt) 3 / 2 (2tt) 3 / 2 



Sp-r(^l) Sr(^2) 



ki.k 2 



3 i(k!+k 2 ).q 



r=p,s=p—r 



4 £ 



p! 



6 ^ (p — r — s)\r\s\ 

r=0,s=0 ' 



r d J ki d J 

7 (2tt) 3 / 2 (27 



3 k 2 d 3 k 3 



)3/2 ( 2yr )3/2 



->p— r— s 



(/Ci)s r (/C 2 )s s (fc3) 



f | ki.k 2 ^ 2 _ ^ k 2 .k 3 y _ ^ k 3 .ki ^ 2 2 ki.k 2 k 2 .k 3 k 3 .ki 



x e i(k!+k 2 +k 3 ).q 



&2&3 



/c 3 /ci 



h 2 h 2 h 2 



(C.20) 

The average of (C.20) over the volume Vb introduces the parameters j p , <i p and the 
functions Ji and J 2 defined in (C.13) and (C.14): 

Jp = ^p + 2^ r! (/-r)! / (27r)3/2 (27r)3 2 /2 ^-r(fci)^(fc2)Ji(k!,k 2 ) 



r=0 
r=p,s=p—r 



+ 



- y 



p! 



d 3 k! d 3 k 2 d 3 k 3 



6 r= f^ =0 (p-r-s)!r!s! y (2vr) 3 / 2 (2tt) 3 / 2 (2tt) 3 / 2 
x J 2 (ki,k 2 ,k 3 ). 



s p _ r _ s (/ci)s r (/c 2 )s s (/c 3 ) 



(C.21) 

which can be integrated out by the use of the properties (C.15) and (C.16). We 
eventually get 



1 r ~ ^ pi 1 ' — 1 

r=0 Vr y 



r=p,s=p—r 



P 



o- / / I 1/ M r s dp-r-s- 

r=0,s=0 ; 



(C.22) 

Then, using the definitions (20) and (C.4) we can derive a relation between Qy and 



g v (t, y)/V = 1 + Q D (t, y) + \g 2 D (t, y) + ^G 3 D (t, y) = [l + Gd (t, y)m s . (C.23) 



The equation (C.3) leads to another relationship of the same kind. The only 
qualitative change is the presence of the operator T(.). But as it is a linear operator, 
it commutes with any operations such as the Fourier transform, ensemble average, 
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and integration over Vq. Similar calculations then lead to 

r=p , 1 r=p,s=p—r 

El T( A \ A ,, 

*>p—r — s 



r=0 v 7 r=0,s=0 

which gives 

T y)\ [1 + y)/3] 2 = -4nGp [V - Q v (t, y)] . (C.25) 



(C.24) 



The equations (C.23) and (C.25) with the definition of T(.) lead to the differential 
equation 

d ' (j? \r u ,A/i/iiM -(Gp + A/47r) Vfr 

since the co moving radius Ro(t) is proportional to the expansion factor and 
satisfies Rq = —4tv/3 J)GRq + ARq/3. The equation (C.26) is the differential 

1 /3 

equation describing a spherical collapse, where Rq [Gv(t, y)/Vo] is the size of the 
object. The behavior of Qv{y) for small values of y ensures that the corresponding 
overdensity is simply y. The function Qv(t, y) is then the volume occupied at the 
time t by an object of linear initial overdensity y and of initial comoving volume 

v . 

This result is not so surprising since we saw that the application of the top- 
hat window function just replaces the geometrical dependences that enter in the 
real dynamics by their means, i.e. , their monopole part, and the dynamics of the 
spherical collapse just corresponds to the case where there is only one monopole. 

C.3 The functions Qyi{t 1 y) and Qs{t, y) 

The purpose of this part is to compute the values of the functions, 

G (t,A 

Gv<t,y) = ^ -a;; -, 

P= ° (C27) 

0s(t>v) = 2^ -2- P Zf. 

p=Q °L V- 

in the rare peak approximation. The first function involves a cumulant between 
p + 2 variables that requires at least p + 1 links between them. It can only be done 
with the moment of 2p + 2 Gaussian variables. As a result 

r=p+l 

<^ 2 (*)^> c « E {V {r \t)V^-r){t)5i) c . (C.28) 

r=l 
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The calculation of such quantities is simple but requires a good understanding of 
the result obtained in the previous part. The order p of the volume can be written 
in such a way, 



(C.29) 



3 3 

viP){t) = V °J • • • (2^7^ • • • 6 K W TR Ilki + • • • + K\ Ro] 

x J(ki,...,kp) 

where J(ki, . . . , k p ) is a function of the angular parts of the wave vectors ki, . . . , k p 
only. The result obtained in the part (C.l) is then a property of the function 
J(ki, . . . , k p ) that reads 



J dfii . . . dn p W TR [|ki + . . . + kp] i? ] J(ki, . . . , k p ) 



( 4yr )p M*) ^ TH (A; lJ Ro)...W TH (A;p J Ro). 



(C.30) 



Indeed the function J(ki, . . . , k p ) turns out to be a combination of the functions J\ 
and J<i for which such a factorization has been explicitly established in the appendix 
B. The relation (C.30) can also be used in (C.28) since the angular parts of the wave 
vectors appearing in the expressions of and of V^ +2 ~^ are independent. As 
a result we simply obtain 



r=p+l 



pi 



" x (r-l)!(p+l-r)! 



jr(t) jp+2-r(t) 



so that 



Gv*(t,y) = ol 



-7-Gv{t,y) 

dy 



(C.31) 



(C.32) 



The same geometrical arguments can be invoked for the calculation of Gi/v(t,y) 
that lead to the expression, 



Gs(t,y) = V /g v (t,y)-l. 



(C.33) 
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APPENDIX D: Perturbative calculations in Lagrangian coordinates 

This Appendix is devoted to the derivation of the third order of the Jacobian 
in a Lagrangian calculation (see the main text for details), and some expectation 
values of interest. The principle of these calculations is not new and has been 
presented elsewhere (see for instance Moutarde et al. 1991, Buchert 1992) 

D.l The three first orders of the displacement field 

To each comoving position q is associated a displacement vector \I/(q) giving 
the present position x of the particle being at q initially, 



x = q + *(q), 



(D.l) 



and the displacement field is thought, in a perturbative approach, to be a small 
correction to the Hubble expansion (contained in q). 

The Jacobian of the transformation between q and x, |<9x/<9q|, gives a 
direct information on the volume occupied by the particles at any stage of the 
dynamics (before shell crossing). Its expression is obtained by the calculation of 
the determinant of the matrix the elements of which are 



= Oij + 



<9q 



dqj 



(D.2) 



where is the Kronecker symbol. It leads to the relation, 



J(t,q) =1 + V q .* + 



(Vq-*) 2 -£**A-.* 



+ 



(V q -*) 3 - 3V q .¥ £ vE^ tf . , + 2 £ * 4J * iifc * fc>i 

ijk 



(D.3) 



where Vq is the gradient taken relatively to q and the summations are made over the 
three spatial components. The divergence of the acceleration field is proportional 
to the density field so that 



J(t, q) ( V x .x + 2- V x -x ) = -AnGp [1 - J(t, q)] (DA) 
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in which a dot denotes a time derivative. In term of the displacement field it reads 
T (V q .*) + Vq T (V q - £ (% ,) + ± (V q .¥) 2 T (V q .*) 

-±T (Vq ]T ^./I',.; - Vq T (* 4J ) J + £ tf^fcT (¥ fc>i ) 



'J 



-47rGp[l- J(*,q)] 



where T(.) is a time derivative operator, 



(0.5) 



T(A) = A + 2-A 
a 



(D.6) 



The displacement field is completely determined by the the non-rotational 
constraint in the x coordinates, which leads to a third equation, 



1 

+ 2 



(Vq.W) 2 -^*,,^,, 



(0.7) 



In the following I will assume that the cosmological parameters are the ones 
of an Einstein-de Sitter universe. 

The principle of the calculation is the following. We assume that the 
displacement field is a small quantity. At a given order the elimination of J(t, q) in 
(D.3) and (D.5) provides an equation for the divergence of the displacement field. 
The equation (D.7) then gives the complete information over its non-symmetric 
part. 

The time dependence of the first order solution, D{t), is solution of the 
linearized equations in the displacement field. As a result D(t) is solution of 



T (£>(*)) = 4nGpD(t). 



(D.8) 



There are two solutions for the previous equations: One is an increasing function 
of time and is proportional to a(t) the other is a decreasing function of time, 
proportional to t . Since we assume that the fluctuations begin to grow at an 
arbitrarily small time the growing mode completely dominates the dynamics. 
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The first order solution then reads 

jU(t, q) = V q q) = £>(f)e+(q), (D.9) 

where e+(q) is a small perturbation in the displacement field corresponding to the 
growing mode. This quantity is related to the initial overdensity that gives birth 
to the cosmic structures. The equation (D.7) insures that \E f | 1 ' ) 
displacement field is completely determined by its divergence 

At the second order the equation for the divergence is 



so that the 



T (Vq .tf (2) ) + Vq .^ (1) (q, t)T (v q 

(v q .*<») 2 -£(*S) ; 



(i) 



47rGp 



13 



(D.10) 



The time dependence of the solution of this equation is D 2 (t) (once the only pure 
growing mode has been selected) . Note that this result is specific to the Einstein-de 
Sitter universe and is not true in general. The exact result for any value of O is 
given by Bouchet et al. (1992). As a result we obtained, 



V q .^ 2 )(t,q) = — 



J (2) M) 



v q .*<") 2 -5:(*s) 2 



l 3 



{D.ll) 



Once again the equation (D.7) insures that ^\ 2 j = ^f}- 

The third order is calculated the same way. Its time dependence is simply 
D(t) 3 and the solution reads 



Vq.^)(t,q) =-± 



9 

-1 
T8 



ij 

(v q .*<») 3 - 3v q .»« y: (»&>)'+ e *&'*}X] 



1 

+ 9 



(v q ,*w) 3 - sv q .»w yj (*w) 2 + e »&'»53»a 

(U."12) 



(3) (3) (3) 

Note that W 4 * V so that a complete determination of ■ cannot be derived 
from (D.12). 



D.2 Calculation of(V@\t)6 L ) 

The second part of this appendix is devoted to the calculation of (V^(t) 8l), 
where V^\t) is the integral of over the volume Vo- In the wave vectors space 
these quantities are given by integrations over k with the filter function Rq) 

as defined in the appendix B. The Fourier transforms, of the first order of the 
Jacobian are defined by 



J«(q) 



d 3 k 



exp(iq.k). 



(2tt) 3 / 2 

As a result the integral of J^ 1 ' over the volume is given by 



(D.13) 



(D.U) 



and 



d 3 k! d 3 k 2 d 3 k 3 

(2tt) 3 / 2 (2tt) 3 / 2 (2tt) 3 / 2 



Xj kJk 2 J k 3 




W TH (|ki +k 2 + k 3 | R ) 



(k!.(k 2 + k 3 )r 



+ 



fc 2 |k 2 + k 3 | 2 

l 2 




k?.k 3 



k 2 .k 3 



+ 2 



k 2 h 

kj.k 2 k 2 .k 3 k 3 .ki 

h 2 h 2 h 2 



+ 



(D.15) 

If the initial density field has Gaussian fluctuations as described in the part 2 the 
random variables are Gaussian and 



7 (1) 
J k 



D(t) 



since S(t, q) = [1 — J(t, q)/J(t, q)] —J^\t, q) at the first order. 

The ensemble average of the product V^'Sl exhibits the expectation value 
of the product of four Gaussian random variables, 



I x •(!)•(!) -(1)\ _ 




(VkAAs) 



^k 6 k lt 



eye. 



(D.16) 



[P(k)5 3 (k + ki) P(A; 2 )5 3 (k 2 + k 3 ) + eye] 
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where eye. stands for three other terms obtained by a circular permutation of the 
indices. The last expression holds because the random Gaussian variables describe 
the density field of a randomly homogeneous and isotropic universe with the power 
spectrum P(k). 

The use of the equations (D. 14-16) gives 

4 



<W 3) 6 L ) = 



Wo 
21 



D(t) 
D{U) 



d 3 k 



(27T) 



x 



1 - 



A^lk + k'l 2 



1 - 



k.k' 

kk' 



(D.17) 



This integral can be simplified since the integration over the angular part between 
k and k' can be made separately. We can define the function /(r) by 



f(r) = 



1 



du _ 

1 1 + t 2 + 2tu 



(l-u 2 f 



(D.18) 



then 



(V^ 6 L ) 



Wo 
21 



D{U) 



d 3 k 



(27T) 



zW^ R (k Ro)P(k) 



3i J 



d^k 



(27T) 



3 P(k')f(k'/k). 



The calculation of /(r) can be done analytically and the result reads 



^ = Tr 



_2A 3 + -A + (1-A 2 ) 2 ln4 ± ^ 
3 v A - 1 



A T 1 

A= 2 + 27- 



(D.19) 



(D.20) 



The shape of the function /(r) is given in Fig. 3. Note that /(r) ~ 5/8 for 
r — > and /(r) ~ 5/8r 2 for r — > oo, so that the integral (D.19) converges only if 
P(k) <C k~ x when k — > oo. 



D.3 Calculation of([V^(t)] 2 ) 

The expression of V^ 2 -*^) derives from the results given in (D.ll), 



D(ti) 



3 3 

/ (2^ ^4>TH(^^o)4>TH(^^) (I-. 2 ) 

(D.21) 



where u = ki.k 2 //ci/c2- As a result we have 



2 r 



D(t) 



xP(h) P{k 2 ) (1-u 2 ) 2 . 



f d 3 kx d 3 k 2 Tir2 ,„ , , „ v 



22) 
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This expression can be simplified by the change of variable 



that leads to 

([v^(t)] 2 ) 



4 
7 

x 



D(t) 
D{U) 



k + = ki + k 2 ; 
k_ = ki - k 2 ; 



d 3 k + d 3 k_ 



(D.23) 



(2tt) 3 (2tt) 



3 Wf H (fc+i2o) P(fci) P(fc 2 ) 



(1 — v 2 ) 



(k+/k-+k-/k+y -4v 2 

(D.24) 

where v = k_|_.k_//c + /c_. In case of a power-law spectrum with an index n (eq. 
[37]) the result reads, 



x 



/>/_ 



g(g) 
+ °° dr r 2+ 



d 3 k + 



k 5 + +2n W 2 (k+R ) 



(27T) 



, (l-t> 2 ) 2 (£;+/£;_ + /c_//e + ) 2 -4?; 2 



n/2-2 

(D.25) 



which has to be integrated numerically. 
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FIGURE CAPTIONS 



Fig. 1 : Evolution of the density profile for an Einstein-de Sitter universe. 
The initial profile is the expected shape of the profile for various value of the 
power-law index: n = — 2 for (a), n = — 1 for (b) and n = for (c). The solid 
lines correspond to the initial profiles. The n = case corresponds to a Poisson 
distribution and the correlation length is which explains this top-hat shape of 
the initial profile. The dashed lines correspond to the exact nonlinear evolution of 
these profiles for various time steps, 5l = 0.25, 0.5, 0.80, 1.10. 

Fig. 2 : The nonlinear density as a function of 5l (eq. [20]) for various 
cosmological models. The solid line is the analytic expression S = (1— 5^/1. 5) -1 - 5 — 1 
corresponding to O — > 0, A = 0; the dotted line is for O = 1, A = 0; and the dashed 
line is for O = 1, A = A/3H$ = 1. 

Fig. 3 : The function /(r) (eq. [D.20]) as a function of r. 

Fig. 4 '■ The corrective coefficient, C e xp (eq. [34]) as a function of n for a 
power-law spectrum (a) and as a function of the cutoff, k c , for a truncated power- 
law spectrum and for various values of n (b). R is the size of the perturbation. 

Fig. 5 : Evolution of the density of a fluctuation as a function of the initial 
overdensity and of the time. The upper solid line corresponds to the spherical model, 
the second solid line corresponds to the spherical model solution truncated at the 
third order. The long dashed and small dashed lines correspond to the expected 
evolution of a fluctuation up to the third order characterized respectively by v = 4 
and v = 2. The corrective coefficients are calculated for n = — 2. 

Fig. 6 : The corrective coefficient as a function of scale in case of a CDM 
spectrum. For the mass scale of the clusters the correction is still small, so that the 
spherical collapse model is expected to be a good description for the collapse of a 
cluster. 

Fig. 7 : The coefficient C sca ^ (eq. [46]) as a function of n (solid line) and 
the coefficient Cg Cat (eq. [48]) (dashed line). 

Fig. 8 : Schematic representation of the situation obtained in the n « — 2 
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case. The expectation behavior (dashed line) is close to the spherical collapse 
model (solid line) with possible fluctuations represented by the shaded area, the 
departure from the spherical collapse and the fluctuations being of the same order 
of magnitude. 



